<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Minimize order of a lowpass FIR filter (magnitude design)</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/filter_design/html/fir_mag_design_lowpass_min_order.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Minimize order of a lowpass FIR filter (magnitude design)</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "FIR Filter Design via Spectral Factorization and Convex Optimization"</span>
<span class="comment">% by S.-P. Wu, S. Boyd, and L. Vandenberghe</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs an FIR lowpass filter using spectral factorization method where we:</span>
<span class="comment">% - minimize the filter order</span>
<span class="comment">% - have a constraint on the maximum passband ripple</span>
<span class="comment">% - have a constraint on the maximum stopband attenuation</span>
<span class="comment">%</span>
<span class="comment">%   minimize   filter order n</span>
<span class="comment">%       s.t.   1/delta &lt;= |H(w)| &lt;= delta   for w in the passband</span>
<span class="comment">%              |H(w)| &lt;= atten_level        for w in the stopband</span>
<span class="comment">%</span>
<span class="comment">% We change variables via spectral factorization method and get:</span>
<span class="comment">%</span>
<span class="comment">%   minimize   filter order n</span>
<span class="comment">%       s.t.   (1/delta)^2 &lt;= R(w) &lt;= delta^2  for w in the passband</span>
<span class="comment">%              R(w) &lt;= atten_level^2           for w in the stopband</span>
<span class="comment">%              R(w) &gt;= 0                       for all w</span>
<span class="comment">%</span>
<span class="comment">% where R(w) is the squared magnited of the frequency response</span>
<span class="comment">% (and the Fourier transform of the autocorrelation coefficients r).</span>
<span class="comment">% Variables are coeffients r. delta is the allowed passband ripple</span>
<span class="comment">% and atten_level is the max allowed level in the stopband.</span>
<span class="comment">%</span>
<span class="comment">% This is a quasiconvex problem and can be solved using a bisection.</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%*********************************************************************</span>
<span class="comment">% user's filter specs (for a low-pass filter example)</span>
<span class="comment">%*********************************************************************</span>
<span class="comment">% filter order that is used to start the bisection (has to be feasible)</span>
max_order = 20;

wpass = 0.12*pi;        <span class="comment">% passband cutoff freq (in radians)</span>
wstop = 0.24*pi;        <span class="comment">% stopband start freq (in radians)</span>
delta = 1;              <span class="comment">% max (+/-) passband ripple in dB</span>
atten = -30;      <span class="comment">% stopband attenuation level in dB</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% create optimization parameters</span>
<span class="comment">%********************************************************************</span>
m = 15*(max_order);   <span class="comment">% freq samples (rule-of-thumb)</span>
w = linspace(0,pi,m);

<span class="comment">%*********************************************************************</span>
<span class="comment">% use bisection algorithm to solve the problem</span>
<span class="comment">%*********************************************************************</span>

n_bot = 1;
n_top = max_order;
n_best = Inf;

<span class="keyword">while</span>( n_top - n_bot &gt; 1)
    <span class="comment">% try to find a feasible design for given specs</span>
    n_cur = ceil( (n_top + n_bot)/2 );

    <span class="comment">% create optimization matrices</span>
    <span class="comment">% A is the matrix used to compute the power spectrum</span>
    <span class="comment">% A(w,:) = [1 2*cos(w) 2*cos(2*w) ... 2*cos(n*w)]</span>
    A = [ones(m,1) 2*cos(kron(w',[1:n_cur-1]))];

    <span class="comment">% passband 0 &lt;= w &lt;= w_pass</span>
    ind = find((0 &lt;= w) &amp; (w &lt;= wpass));    <span class="comment">% passband</span>
    Ap  = A(ind,:);

    <span class="comment">% transition band is not constrained (w_pass &lt;= w &lt;= w_stop)</span>

    <span class="comment">% stopband (w_stop &lt;= w)</span>
    ind = find((wstop &lt;= w) &amp; (w &lt;= pi));   <span class="comment">% stopband</span>
    As  = A(ind,:);

    <span class="comment">% This is the feasiblity problem:</span>
    <span class="comment">% cvx_begin quiet</span>
    <span class="comment">%     variable r_cur(n_cur+1,1);</span>
    <span class="comment">%     10^(-delta/10) &lt;= Ap * r_cur &lt;=  10^(+delta/10);</span>
    <span class="comment">%     abs( As * r_cur ) &lt;= +10^(+atten/10);</span>
    <span class="comment">%     A * r &gt;= 0;</span>
    <span class="comment">% cvx_end</span>
    <span class="comment">% Unfortunately it seems to be a bit unreliable to solve. So we have</span>
    <span class="comment">% reformulated it as a stopband minimization. If the optimum stopband</span>
    <span class="comment">% attenuation is smaller than 10^(atten/10), then we have feasibility.</span>
    <span class="comment">% formulate and solve the feasibility linear-phase lp filter design</span>
    cvx_begin <span class="string">quiet</span>
        variables <span class="string">r_cur(n_cur,1)</span>;
        minimize( max( abs( As * r_cur ) ) );
        10^(-delta/10) &lt;= Ap * r_cur &lt;= 10^(+delta/10);
        A * r_cur &gt;= 0;
    cvx_end

    <span class="comment">% bisection</span>
    <span class="keyword">if</span> isnan( cvx_optval ),
        fprintf( 1, <span class="string">'Solver failed for n = %d taps, assuming infeasible\n'</span>, n_cur );
        n_bot = n_cur;
    <span class="keyword">elseif</span> cvx_optval &lt;= 10^(atten/10), <span class="comment">% strfind(cvx_status,'Solved') % feasible</span>
        fprintf(1,<span class="string">'Problem is feasible for n = %d taps\n'</span>,n_cur);
        n_top = n_cur;
        <span class="keyword">if</span> n_best &gt; n_cur, r = r_cur; <span class="keyword">end</span>
    <span class="keyword">else</span> <span class="comment">% not feasible</span>
        fprintf(1,<span class="string">'Problem not feasible for n = %d taps\n'</span>,n_cur);
        n_bot = n_cur;
    <span class="keyword">end</span>
<span class="keyword">end</span>

h = spectral_fact(r);
n = length(r);
fprintf(1,<span class="string">'\nOptimum number of filter taps for given specs is %d.\n'</span>,n);


<span class="comment">%********************************************************************</span>
<span class="comment">% plots</span>
<span class="comment">%********************************************************************</span>
figure(1)
<span class="comment">% FIR impulse response</span>
plot([0:n-1],h',<span class="string">'o'</span>,[0:n-1],h',<span class="string">'b:'</span>)
xlabel(<span class="string">'t'</span>), ylabel(<span class="string">'h(t)'</span>)

figure(2)
<span class="comment">% frequency response</span>
H = exp(-j*kron(w',[0:n-1]))*h;
<span class="comment">% magnitude</span>
subplot(2,1,1)
plot(w,20*log10(abs(H)),<span class="keyword">...</span>
    [wstop pi],[atten atten],<span class="string">'r--'</span>,<span class="keyword">...</span>
    [0 wpass],[delta delta],<span class="string">'r--'</span>,<span class="keyword">...</span>
    [0 wpass],[-delta -delta],<span class="string">'r--'</span>);
axis([0,pi,-40,10])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'mag H(w) in dB'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
Problem not feasible for n = 11 taps
Problem not feasible for n = 16 taps
Problem is feasible for n = 18 taps
Problem is feasible for n = 17 taps

Optimum number of filter taps for given specs is 17.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_mag_design_lowpass_min_order__01.png" alt=""> <img src="fir_mag_design_lowpass_min_order__02.png" alt=""> 
</div>
</div>
</body>
</html>